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ABSTRACT 

We present a new set of 3.5 Post-Newtonian equations in which Newtonian hydrodynamics is coupled 
to the nonconservative effects of gravitational radiation emission. Our formalism differs in two signifi- 
cant ways from a similar 3.5 Post-Newtonian approach proposed by Blanchet (1993, 1997). Firstly we 
concentrate only on the radiation-reaction effects produced by a time-varying mass-current quadrupole 
Sij. Secondly, we adopt a gauge in which the radiation-reaction force densities depend on the fourth 
time derivative of Sij, rather than on the fifth, as in Blanchet's approach. This difference makes our for- 
malism particularly well-suited to numerical implementation and could prove useful in performing fully 
numerical simulations of the recently discovered r-mode instability for rotating neutron stars subject to 
axial perturbations. 

Subject headings: relativity: Post-Newtonian approximation — gravitational radiation reaction — 
stars: neutron — stars: oscillations — instabilities. 

1. INTRODUCTION 

The recent discovery of an instability in r-mode oscillations in relativistic rotating stars has generated widespread 
interest, r-mode oscillations in Newtonian rotating stars have been widely investigated in the past (see, for instance, 
Papaloizou and Pringle 1978; Provost et al. 1981; Saio 1982), but the evidence that they are indeed unstable to the 
emission of gravitational radiation is rather recent. The first numerical calculations carried out by Anderson (1998) and 
confirmed analytically by Friedman and Morsink (1998) have spawned a growing literature on the subject (Andersson, 
Kokkotas and Schutz 1998; Kojima 1998; Kokkotas and Stergioulas 1998; Levin 1998; Lindblom and Ipser 1999; Lindblom, 
Mendell and Owen 1999; Lindblom, Owen and Morsink 1998; Lockitch and Friedman 1998; Madsen 1998; Owen et al. 
1998). Much of the interest in r-mode oscillations is related to the fact that their existence is not dependent on a specific 
rate of rotation; these modes are, in fact, unstable for arbitrarily slowly-rotating, perfect fluid stars. This result represents a 
significant difference from previously investigated, rotation-induced instabilities, as for instance the bar-mode instability, 
which requires a minimum rotation rate of the star (Chandrasekhar 1970; Friedman and Schutz 1978; Lindblom and 
Mendell 1995; Stergioulas and Friedman 1998). Consequently, the r-mode instability may have a more pervasive effect. 

The r-mode instability is a purely relativistic effect, triggered by the emission of gravitational radiation and can be 
explained in terms of the basic Chandrasekhar-Friedman-Schutz instability mechanism (Chandrasekhar 1970; Friedman 
and Schutz 1978). Since gravitational radiation removes positive angular momentum from a prograde mode (i.e. a mode 
that in an inertial frame is seen as moving in the same positive Lp direction as the star) , it will also extract positive angular 
momentum from any perturbation which (as a result of the star's rotation) is prograde in the inertial, but retrograde 
in the corotating frame. Such a mode has, in the corotating frame, negative angular momentum (the perturbed fluid 
does not rotate as fast as it did without the perturbation) and by making its angular momentum increasingly negative, 
gravitational radiation drives an instability (Friedman 1998). A significant difference between the axial r-mode instability 
and the previously known gravitational radiation driven bar-mode instability, is that gravitational radiation couples with 
r-mode oscillations primarily through time- varying mass-current multipole moments rather than through the usual time- 
varying mass multipole moments. We here present a new set of Post-Newtonian (PN) hydrodynamical equations which 
make the calculation of radiation-reaction forces due to mass-current multipole moments numerically feasible. 

There are a number of reasons why a numerical investigation of the onset, growth and saturation of the r-mode 
oscillations is of great interest. All investigations to date have been based on perturbation analyses usually truncated at 
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the lowest order in the expansion parameter, namely the star's angular velocity (Andersson 1998; Friedman and Morsink 
1998; but see Andersson, Kokkotas and Schutz 1998; Kokkotas and Stergioulas 1998; Lindblom, Mendell and Owen 1999 
for a treatment at higher order). Within this framework, the perturbation is parameterized in terms of the angular velocity 
of the star and its growth is heuristically followed in terms of the energy and angular momentum losses produced by the 
gravitational radiation emission (Owen et al. 1998). The relevant timescales for the growth and the subsequent viscous 
decays of the instability are also estimated through perturbation analysis (Lindblom, Owen and Morsink 1998). Although 
neglecting the presence of a magnetic field, these approaches have clarified the basic features of the instability and provided 
the first estimates of the importance of the instability in extracting angular momentum from hot young neutron stars, thus 
setting an upper limit on their angular velocity (Lindblom, Owen and Morsink 1998; Andersson, Kokkotas and Schutz 
1998). They also provided qualitative and quantitative information about the expected gravitational waveforms (Owen 
et al. 1998). However, they are not capable of describing the nonlinear development of the instability or identifying its 
point of saturation. For these features, numerical simulations are required. 

To study this mode numerically, it will be useful and adequate to follow the conservative Eulerian hydrodynamics via 
Newtonian equations and all the nonconservative effects due to (mass-current multipoles) gravitational radiation emission 
via a Post-Newtonian radiation-reaction potential at the 3.5 order [this is what we shall refer to as (0 -|- 3.5) PN]. Such a 
(0-1-3.5) PN approach has the advantage of capturing of all the relevant nonconservative general relativistic effects without 
having to resort to a more complicated relativistic hydrodynamics treatment. Moreover, a (0-1-3.5) PN treatment allows us 
to clearly disentangle the different sources of gravitational radiation. As a result, in a simulation of the r-mode instability, 
we can selectively neglect all the dissipative contributions coming from mass multipole moments and concentrate solely on 
the mass-current quadrupole moment, which we expect to be the dominant mass-current multipole moment. Disentangling 
these modes is not possible in a full general relativistic treatment. Finally, because simulations of the onset and growth 
of the r-mode instability also require the numerical evaluation of stable configurations on growth timescales much longer 
than the dynamical timescales (i.e. growth timescales ^ rotation period), we also expect that a three-dimensional, fully 
relativistic simulation may be, at present, beyond reach. On the other hand, a radiation-reaction formalism, though 
approximate, allows us to use a scaling in order to artificially accelerate the onset of an instability without changing the 
underlying physical evolution. 

The organization of the paper is as follows: in Section ^ we summarize the basic steps of a PN expansion and the 
hydrodynamical equations that emerge from it. Particular attention will be paid to the radiation-reaction forces and 
to the losses they produce in the energy and angular momentum of the system. In the following Sections ^ and ^ we 
discuss two explicit expressions for the radiation-reaction force due to time-varying mass current quadrupole moments. 
The first one was obtained by Blanchet (1997) (Section ||), while the second is derived in a new gauge (Section ^ which 
makes it better suited for numerical implementation. Section ^ also contains detailed verifications that the new force 
yields the required rates of energy and angular momentum loss. Section ^ synthesizes the main results derived in the 
previous Sections and, for the benefit of the reader interested in numerical implementation, presents the final form of the 
hydrodynamic equations and radiation-reaction terms. Having in mind the study of the r-mode instability, we present 
in Section ^ a useful rescaling of the radiation-reaction term which will accelerate the growth-time of the instability to a 
timescale set by numerical constraints. Section ^ summarizes our conclusions and the prospects for numerical computations 
exploiting the formalism introduced here. Throughout the paper we will adopt Cartesian coordinates. G and c denote 
the gravitational constant and the speed of light. Greek indices run from to 3, Latin indices from 1 to 3, and we use 
Einstein summation convention on matched indices. 



2. PN EXPANSION: THE BASIC EQUATIONS 

In this Section we briefly summarize the standard approach to perform a PN expansion for the equations of relativistic 
hydrodynamics. (A recent discussion of the foundations and applications of the PN approximation has been given by 
Asada and Futamase 1997). We adopt the standard 3-1-1 splitting of spacetime and write the line element in the form 

ds^ = g^^tlxf^dx" = -{a^ - P^f3')c^dt^ + 2(3,cdtdx' + -/tjdx'dx^ , (1) 

where a, /3% and jij are the lapse function, the shift vector, and spatial 3-metric, respectively, while f3i = jijP-' ■ 
We also consider a perfect fluid, whose energy-momentum tensor is 

T^'^ = {pc^ +pe + P)u''u'' + Pg'''' , (2) 

where p is the rest-mass density, e is specific internal energy, P the pressure, and u'^ the fluid four-velocity. It is also 
convenient to introduce a coordinate velocity u*, defined as 

= + (3) 



Imposing the conservation of rest-mass, energy and momentum yields the standard relativistic hydrodynamic equations. 
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(7a) 
(7b) 



and 7 = det(7ij). Equations (§)"(|3) will also be referred to as the continuity, energy and Euler equations, respectively. 

We next proceed to the PN approximation and perform a series expansion of the metric in the inverse powers of c up 
to the 3.5 PN order (Chandrasekhar 1965; Chandrasekhar and Esposito 1969; Asada, Shibata and Futamase 1996): 



1 i 1 1 1 1 1 1 

a = 1 + —(j) + —iCt + —^a + —&a + —ja + —9,a + — ga + 0(c 

c° c° c' 



(8a) 
(8b) 
(8c) 



where the left subscript n indicates the coefficient of the 0(c~") term in the series expansion, and (p is the Newtonian 
potential. We implicitly adopt the usual PN gauge in which 20. = and 2hij — —20(5^; (Chandrasekhar 1965; Chan- 
drasekhar and Esposito 1969) and 5a is a function of time only. Similarly, the four-velocity is also expanded in terms of 
u*, and full expressions for this are given by Chandrasekhar (1965), Chandrasekhar and Esposito (1969), Asada, Shibata, 
and Futamase (1996). 

Using (^), the 3.5 PN expression of the Euler equation (^ can then be written as 
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Hereafter we neglect the higher order difference between contravariant and covariant components and use the latter only. 

It is convenient to distinguish, in the right hand side of (|^), the terms related to dissipative radiation-reaction effects 
from those arising from conservative hydrodynamical stresses. In particular, we define the radiation-reaction force densit}^ 



prr _ p2.5PN 



:,3.5PN 



rrm ^ prrc ^ (^q^ 

where F''"" and F''"^ refer to the radiation-reaction force due to time-varying mass multipole moments and mass-current 
multipole moments, respectively. The general (slice-independent) expressions for ir-S-SPN ^3.5PN given by (Asada 
and Futamase 1997; Blanchet 1997) 



^-ij^2.5PN ^ _ ^ Vjd^iePj) - v,d,{ef3,) 

-dtishijVj) - zhijVkdkVj (11) 

^-lj^3.5PN ^ _ _ + dd,8l3j)vj 

-dtijhijVj) - VkdkirhijVj) + ^VjVkdiirhjk) + SPf-'^^^ira, ePi, bhij) , (12) 

where di — d/dxi^ dt = d/dt, and in deriving (^ij) we have exploited that 9i[5Q;(t)] = 0. Note that while F^-^^^ 
is dependent only on the time-varying mass quadrupole moments , F^'^^'^ is, in general, dependent on time-varying 
mass quadrupole moments, mass octupole moments, as well as on time-varying mass-current quadrupole moments. In 

^Hereafter we refer to the radiation-reaction force densities simply as radiation-reaction forces. 
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particular, in equation ([l^), we have symbolically indicated by (5F'^'^^^(7a, g/3i, ^hijVj) all of the contributions coming 
from mass quadrupole moments. 

The work done and the torque produced by the radiation-reaction forces F''' must balance the energy and angular 
momentum carried off to infinity by the gravitational waves. The emission rate of gravitational waves is known from the 
multipolar decomposition of the radiation field (Thorne 1980). In the absence of any dissipative mechanism other that 
gravitational wave emission, the energy and angular momentum loss rates can be readily calculated as 

(13) 

^= js^e,,,x,Fl\ (14) 

where d/ dt ~ dt + vidi and e^-fc is the Levi-Civita symbol. The total energy and angular momentum have the usual 
Newtonian definitions. 




with = VkVk- Equations ( |T^ ) and (Q) will be used repeatedly in this paper to verify the correctness of the derived 
expressions for the radiation-reaction forces. 

So far, our discussion of radiation-reaction forces has been general and we have not restricted ourselves to a specific 
physical configuration and considered a generic gauge in which the expansion (H) is valid. Hereafter, however, we will 
concentrate on a PN formulation of the radiation-reaction forces which may prove useful in a numerical study of the 
r-mode instability. Because the r-mode instability is predominately excited by mass-current quadrupole moments (cf. 
Section Q), we will neglect any contribution to the radiation-reaction forces coming from mass multipole moments and 
consider F"''° = = (SF'^-^^'^. Even with this restriction, the numerical computation of F'"''^ in presently adopted gauges 
(Burke 1971, Blanchet 1997) is nontrivial. In the following two Sections we discuss these complications and offer a way 
to simplify them. 

3. MASS-CURRENT QUADRUPOLE MOMENT RADIATION-REACTION; BLANCHET'S GAUGE 

The first expression for the radiation-reaction force due to time- varying mass-current quadrupole moments was derived 
by Burke using a matched asympotics expansion (Burke 1969, 1970, 1971) and expressed in terms of vector spherical 
harmonics. His expression, however, does not yield the required energy and angular momentum losses (see Walker and 
Will 1980 for an explanation of the error in the formulation). More recently, a new complete treatment of the radiation- 
reaction and balance equations at 3.5 PN order has been provided by Blanchet (1997) as an extension of earlier work on 
gravitational radiation-reaction forces (Blanchet 1993). Here we briefly review the key steps necessary for our modified 
treatment presented in the next Section. Firstly, the "canonical form" of the linearized metric hf^'^^ is constructed in the 

harmonic gauge (Thorne 1980). This linearized metric expresses the linear deviation from the Minkowski metric 77^" in a 
series expansion in G [i.e. h'^'^ = \/—gg^'^ — rj^^'^ = h'^^-^ + h'^^^ + 0{G^)]. Then, h'^^,^ can be rewritten in terms of two (and 

only two) sets of time- varying multipole moments, referred to as the "mass- type" and "current-type" moments (Thorne 
1980; see Appendix B for details). The contributions to the linearized h'^^^ coming from radiation-reaction effects are 

then derived by taking the half-sum and the half-difference of the retarded and advanced expressions of the multipole 
moments, and by studying the nonlinear corrections by means of a Post-Minkowskian method (Blanchet 1993, 1997). In 
doing this, an infinitesimal gauge transformation to the "generalized Burke-Thorne gauge" (hereafter, we will refer to it 
simply as Blanchet's gauge) is performed in order to obtain h'^'^ , a simplified form of the metric (Blanchet 1993, 1997). 

Restricting our attention only to the radiation-reaction force produced by a time-varying mass-current quadrupole 
moment, it is clear from (|l2|) that we need explicit expressions for the metric coefficients ga, g/^fe and 7/1^ . While the last 
two are known already from the linear term of ft,^", the first one needs to be obtained through an iteration involving also 
the nonlinear terms. As a result of this iteration, the relevant parts of expanded metric functions are [cf. equations (3.6) 
of Blanchet 1997] 



ga = , 

_ 16 G (5) 
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(17a) 
(17b) 
(17c) 
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where Sij is the Newtonian mass-current quadrupolc moment defined as 



(18) 



It is useful to remark that Sij is trace- free, i.e. Sij 6^^ = 0. The right superscript (n) indicates the n-th total time 
derivative: 



(19) 



and = ^{Aij + Aj,). 

Using (|l^), the contribution to the 3.5 PN radiation-reaction force due to a time- varying mass-current quadrupole 
moment in Blanchet's gauge is given by 



-1 jprrc 
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45' 
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(20) 



The validity of expression (p(]|) can be verified by computing the energy and angular momentum dissipation rates. In 
particular, inserting ( |20| ) into (]T^), we immediately obtain 



dE 16 G 
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(21) 



Assuming nearly-periodic motion of the matter field, and averaging over several periods, we can discard the total time 
derivative term and obtain the standard formula of the energy loss due to mass-current quadrupole radiation (Thorne 
1980): 
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(3) q(3)\ 
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where, as usual, 
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Similarly, by using expression (20) in equation (|14D, we obtain 

d^x p {x^XjXkSff} - \yi\^XjS\f) + 2{xjXkViSff} - XjXkVkS'^'') 
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(23) 



(24) 



where |xp = aJia;^. After averaging ( p4[ ) over several periods (and taking two integrations by parts), the formula for the 
angular momentum loss is 



dSi 
'dt 



32G 

'45c 



1..., /^(a) c(3)x 



(25) 



which is again in agreement with the required expression (Thorne 1980). 

Although Blanchet's formalism is clear and complete [indeed Blanchet (1993) and (1997), also discusses 3.5 PN radiation- 
reaction potentials due to mass quadrupole and octupole moments] , it is not particularly simple for numerical implemen- 
tation. The reason for this is evident from expression (pO|), in which the radiation-reaction force depends on a high time 
derivative of the mass-current quadrupole moment Sij. It is often possible, and highly convenient in a numerical calcu- 
lation, to replace some of the time derivatives of the mass and mass-current multipole moments by quadratures. This 
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involves combining the continuity and Euler equations and introducing some supplementary variables (most notably the 
partial time derivatives of the Newtonian gravitational potential) which satisfy elliptic equations (Nakamura and Oohara 

1989). Here, however, obtaining these quadratures is realistic for Slj at most (see Appendix A for a discussion), leaving 

(3) 

the higher order time derivatives to be obtained via finite differencing of Slj . The latter operation can be extremely 
inaccurate, even for a numerical scheme which is second order accurate in time and space, and might introduce numerical 
instabilities. 

To overcome this difficulty, we employ a different gauge condition and derive an alternative form of the metric in the 
next Section. This alternative leads to radiation-reaction forces which are dependent only on the fourth time derivative 
of Sij, a considerable improvement computationally. 

4. MASS-CURRENT QUADRUPOLE RADIATION-REACTION: A NEW GAUGE 

In this Section we adopt a gauge choice different from Blanchet's to obtain a more desirable form for the radiation- 
reaction force. As in Section (^), we start with the linear metric in the canonical form h'^^-^ (Thorne 1980; Blanchet 1993) 

and define h'^^^ = h'^^^ — ^rj^'^hi^i-j, where = (/i(i))'^^. In this gauge, the metric coefficients sf^k and 7/iy are (see 
Appendix B for details) 

4 G 

= = -—e^jkXjXiS[^\ (26a) 
7% = = --^^k<^ki{iS^)i ■ (26b) 

To eliminate the dependence on the fifth time derivative of Stj , we perform an infinitesimal gauge transformation 

hf^i, = hf^i, + d^if, + df,^^ , (27) 
where /i'"'^ = /i'"'^ - \r^f"'h, with h = h''^. We use the freedom in the gauge transformation to set sPk to zero by choosing 



Co 



4G 
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(28a) 
(28b) 



which then yields 



8/3. = , (29a) 

32 G r^f4) /r\r\i \ 

rhij = — ^Xkeki{i^f)i ■ (29b) 

We still have not determined ga, but this can be done by choosing a time-slice condition. Note, however, that selecting 
a specific form for the shift and the spatial 3-metric through equations (|2^) restricts the set of possible choices for 9 a (see 
Appendix C for a discussion). In particular, we impose the maximal slicing condition, for which the trace of the extrinsic 
curvature tensor is set to zero (Arnowitt, Deser and Misner 1962; Smarr and York 1978; Schafer 1983; Blanchet, Damour 
and Schafer 1990) and which is compatible with conditions (p9|). This results in a linear elliptic equation for the lapse 
function, whose 3.5 PN approximation is (see Appendix C for details) 

A(9a) = irh^.d,^) , (30) 
where A denotes the flat spatial Laplacian. Introducing a scalar "superpotential" x satisfying (Chandrasekhar 1969) 

Ax ^2^, (31) 
and using the fact that di^-jhij) = = A(-/hij), we then obtain 

ga = i irh.jdijx) ■ (32) 
The expression for ga can be further simplified by solving ( ^ ) for x 

x(x)^-Gy'dVp(x')|x-x'| . (33) 
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When ( |33[ ) is inserted in equation ([S^), it leads to 



--jG{7hij] 



' dxi 



/ P(x') 9 /"^a^/ 



where rhijSij ~ 0. The vector potential P in (R2|) is defined by 



P(x) eeG J d^x' p(x')| 



(34) 



(35) 



|x - X' I 

and can be most easily calculated by solving the linear elliptic equations 

APi = -AnGpx, . (36) 

While solving equations ( ^6| ) represents an additional computation, nonexistent in Blanchet's formulation, this integration 
is generally not too taxing in a numerical hydrodynamical simulation which already must solve Poisson's equation for the 
Newtonian gravitational potential (Nakamura and Oohara 1989; Oohara and Nakamura 1990, 1991; Shibata, Nakamura 
and Oohara 1992; Ruffert, Janka and Scha fer 1996). 

Finally, using expressions (|4|) and ( |29b| ) for ga and y/iy in (|l2|), we obtain the new gauge expression for the radiation- 
reaction force due to a time- varying mass-current quadrupole moment 



p ^F"" = -di{ga) - dtirhijVj) ~ Vkdkirh^jVj) + ^vjVkdiijhjk) 



(37) 



In subsections 4.1 and 4.2, we verify ( |37| ) by computing the energy and angular momentum loss rates. The reader 
wishing to omit this discussion may proceed directly to Section ^. 

4.1. Rate of Energy Loss 

We calculate the rate at which the total energy of the system is lost to radiation, by substituting expression ( |37| ) into 
equation ([l3|). The relevant integrals that emerge are 



d^x pvtdt{ga) = - / d^x di{pvi)ga = / d^x dt{p)9a 
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d^x pviVjVkdkijhij) = 



(41) 



In derivin g (pa )-(g|) we have made use of the Newtonian continuity and Euler equations as well as of the relation [cf. 
equation (|76a|) in Appendix A] 



Sl^^ = J d^x ptu(iVj)XkVi + JJ d^x d^x' p{x)p{x!) 
Grouping all the terms, we therefore obtain 



^kl(iXj)XkXl 



(42) 
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dE_16Gf (1) (5) £.(2) ^(4) ,.(5) f /•^3„ ^3„/ XjXkx'i 

(fx (fx! p{x)p{x')- 



16 G 



45 c7 ^ 45 c7 dt 



16 G d 



r'|3 



(43) 



As done in Section we now assume quasi-periodicity in the mass-current quadrupole moments and average expression 
43) over several periods. This allows us to discard the total time derivative term and finally obtain the required result 

!)■ 

4.2. Rate 0/ Angular Momentum Loss 
Using equation (p7|), equation (p^ leads to the following terms: 



,3 ^16^ 
d X peijkXjOkga 
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^/|3 il 



^/|3 bk 



i J d^X p{x)eijkXjViV,ndk7hl„i = ^G J d^X p{xbViVmSl,'^ - XlVlV,nS'l2) , 

where we used the Newtonian equation of motion to derive the right-hand side of equation (|4^). 
Grouping all the terms, we then obtain 
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dt 



|x-x\ 
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^'13 li 



(44) 



(45) 
(46) 



(47) 



Once again, we can average equation (|4^) over several periods and obtain the required result (|25D. 

5. (0 + 3.5) PN HYDRODYNAMIC EQUATIONS 

In this Section we present the final set of (0-1-3.5) PN hydrodynamical equations in which the 3.5 PN radiation-reaction 
forces depend only on a time- varying mass-current quadrupole moment. We will present them in a general gauge first 
and then distinguish the expressions resulting from Blanchet's gauge and from our new gauge. 

The general expressions (|) and (|^) for the continuity and Euler equations can be rewritten as 



dp_ _^ djpvi) 
dt dxi 

d{pwk) , d{pwkVi) 



, 



dt 



dxi 



dxk ^ dxk dxk 



2 dxk 



where 



(48) 
(49) 

Wk = cuk = e(8/3fc) + Vj[5jk + e{jhjky\ , (50) 

and where we define e = to highlight the radiation-reaction contributions. Note that the left hand side of equation 
dig] ) contains a partial time derivatives of Wk rather than of Vk- Doing this removes the partial time derivatives of 8/3fe 
and jhjkVj from the right hand side [cf. equation (|T^)]. Moreover, since Wj — Vj = 0(e), all the quantities wj on the right 
hand side of equation (Eoh can be replaced by the equivalent quantities Vj, whenever this is numerically more convenient. 
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In a similar way, the energy equation (||) can be rewritten as 



dt dxi dx. 

or, if we define E = e + hwkWk, in the equivalent form 



d{pe) ^ d{pevi) ^ _pdvi_ ^^^^ 



d{pE) d{pE + P)vi 



dt dxi 

(52) 

In the specific case of an equation of state 

P^{T- l)pe , (53) 
the energy equation (|5^ ) can be written in a (third) simpler form 

where e = {pe)^/'^ . 

As discussed in Sections ^and^ the metric coefficients ga, g/3fe and y/iy appearing in (^)-(p^) represent the radiation- 
reaction potentials and their expressions vary according to the gauge assumed. In particular, they have the form 



Afew Gauge 



1 



sPk = , 



Blanchet' s Gauge 
ga = , 



« 16 G (5) 
sPk = —r^etjkXiXib^ 



rhij — 



32 G 



45 



(4) 



45 



rhrj = , 



(55) 



where Pi is the solution of equation (|36|). We stress that the most important difference in the two gauges is given by the 
appearance of a fourth or of a fifth time derivative of Sij. Note also that, in the new gauge, both the last term of equation 
(|52|), as well as terms including sPk, are zero. 

Finally, the set of hydrodynamical equations is closed by the Poisson's equation for the Newtonian gravitational potential 



A0 = AttGp . (56) 

In the case of the new gauge, equation (|5^) needs to be supplemented by three additional elliptic equations for the 
components of the vector potential Pi 



Boundary conditions at r 



APi = -AnGpx^ . 
oo for the linear elliptic equations (|5l) and (p7|) are given by 



G 



P 



Gxk 



cPyi XkX.,p + 0{r ^) 



(57) 



(58) 



Computing the amplitude and waveforms of the gravitational waves emitted is clearly of great interest since they provide 
the contact with the observations and can be used to extract astrophysical information about the source. In the wave 
zone (Thorne 1980) and at a distance r = |x| from the source (where r ^ L, with L being the size of the source) the 
gravitational wave field is described by the transverse-traceless (TT) part of the linear 3-metric perturbations (Thorne 
1980; Blanchet 1993, 1997). 



{h^^r^it, x) 



4G^ 

^ 1=2 



(-l)'l 



dL-2 



(2) 



ijL-2 



TT 



+ 0(r- 



(59) 



10 



Constructing a Mass-Current Radiation-Reaction Force 



where Ol = diidi2 ■ ■ - di^ and Al^lt — r/c), SL{t — r/c) are the {L — iii2 ■ ■ ■ ii)-th mass muhipole moment and mass-current 
multipole moment respectively. The superscript TT refers to projecting out the transverse, traceless part: 



[A 



iTT 



I Pj ni / m 



}_p p A 

2 ij I'l'n I'm 



(60) 



riirij is the projection operator and 6ij the usual Kronecker-delta symbol. Restricting our attention to 



where = 5. 

the contribution given by mass-current quadrupoles, equation (M) then gives 



c(2) o(2) 



rib 
r 



(61) 



Note that at the 3.5 PN order, there are no gravitational-wave "tail effects" and therefore the mass-current quadrupole 
Sij corresponds to the gravitational wave moment observed. The usual states of polarization of the gravitational waves 



emitted, h+{9, (p) and hx{d, ^p) at a coordinate position (0, ip) on the 2-sphere of radius r, can be derived from (61) 
after selecting the orientation of the source and thus the direction of propagation of the waves (Rasio and Shapiro 1994, 
Ruffert, Janka and Schafer 1996). 

6. TIMESCALES AND RESCALING 

As mentioned in Section we here further explore the possibility of using the set of (0 + 3.5) PN hydrodynamical 
equations presented above to investigate the onset and growth of the r-mode instability. In particular, we want to address 
the problem of the timescales and propose a strategy to suitably rescale them. 

It is commonly assumed that the evolution of the r-mode instability in a unmagnetized, rotating neutron star proceeds 
through three stages (Owen et al. 1998). During the initial stage, any infinitesimal (axial) perturbation grows expo- 
nentially in a timescale t^^^, set by gravitational radiation-reaction. This is followed by the intermediate stage during 
which the amplitude of the mode saturates due to (not yet well understood) nonlinear hydrodynamic effects; the star is 
progressively spun-down as a result of the angular momentum loss via gravitational waves. The Enal stage of the evolution 
occurs when the star's angular velocity is so small that viscous dissipative effects dominate the radiation-reaction forces 
and the r-mode oscillations are damped out. The first stage, for a £ = m = 2 mode, has been estimated to be of the order 
of a few seconds for a neutron star initially rotating at the break-up limit for several different equations of state, while 
the second to be of the order of about one year (Lindblom, Owen and Morsink 1998). 

One complication in simulating r-mode oscillations is that the "natural" timescale r^j^ for the instability to grow and 
saturate is likely to be much longer than the timescale over which a numerical computation can be carried out. Even the 
most sophisticated three-dimensional Newtonian numerical codes suffer from numerical viscosity and are able to preserve 
accurate configurations only for a limited number of stellar rotations (i.e. ^ 10 — 100) and this might well be insufficient 
for the instability to saturate. Below we review the relevant timescales for the r-mode instability and propose a strategy 
whereby, with suitable scaling, we can achieve these timescales in a numerical simulation. Our brief review follows closely 
the results presented by Lindblom, Owen and Morsink (1998). 

Perturbation analysis can be used to estimate t^^ by assuming that the rate of energy loss to gravitational radiation 
emission grows according to 




2E 



(62) 



where £^ > is the energy in the mode and < 0. In the corotating frame of the equilibrium unmagnetized star, E 
can be calculated as 



E 



Sp 



5(j) 6p* 



(63) 



The lowest order expressions for the Eulerian density perturbation 6p and velocity perturbations Sv"" can be deduced 
from the perturbed fluid equations and, in a spherical coordinate system (r, 9, tp), have the form (Lindblom, Owen and 
Morsink 1998; Lindblom, Mendell and Owen 1999) 



5p 



(2£+ 1) fa^R^n'^\ dp 



2£ 



2e+l \ £ + 



/ r \ 



P e{e + i)V2Fr3 V V2i + s J dp 

where (55* (r) is proportional to the gravitational potential, and the axial velocity perturbations are given by 



RJ 



(64) 



(65) 



Here, R and are the radius and angular velocity of the unperturbed star, a^(t) is a dimensionless coefficient parameter- 
izing the amplitude of the perturbation, uj is the (Eulerian) frequency of the mode, and is the magnetic-type vector 
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spherical harmonic. Given an axial perturbation with periodic dependence gH™¥'+'^t) a,nd the definition ( |6^ ) of the energy 

in the mode, the perturbed fluid equations yield the following general expression for the time derivative of E (Ipser and 
Lindblom 1991; Lindblom, Owen and Morsink 1998) 



|2 



(66) 



where r] and C are the shear and bulk viscosities of the fluid (taken as given functions of the density and temperature) 
and 



>(^- 1)[(2^+1)!!]2 

It is easy to distinguish in expression ( |66|) the suppressing viscous terms, driven by the perturbed shear 5aab and 
expansion Sa, from the driving gravitational radiation terms, driven by the time-varying mass 5Ii,m and mass-current 
SSirn multipole moments of the perturbed fluidQ. The explicit contribution to the imaginary part of the frequency of the 
mode due to gravitational radiation-reaction can then be calculated as (Lindblom, Owen and Morsink 1998) 

' I 1 ' pr-^'^+^dr . (68) 




where, as first pointed out by Papaloizou and Pringie (1978), we have used the following relation between the frequency 
of the mode and the star angular velocity 

2mn (^-1X^+2) 
uj = uj.ot -mn^ £{£+!) " "^^ = ^ ■ (69) 



Here Wjot is the angular frequency of the mode in the corotating frame and the last expression in (p9D refers to the 



case in which m = £. Note that the contribution to the growth rate in (68) comes solely from the current multipole 
moments SS£i since we have implicitly neglected the contributions coming from the mass multipole moments Slgg. Such 
an approximation is reasonable because the density perturbations are one order in Q higher than the correspondent velocity 
perturbations and because the density perturbations generate gravitational radiation at a higher frequency (Lindblom, 
Owen and Morsink 1998). 

The general expression for r^j, in (|6^) can be rewritten in a number of alternative ways, some of which are more useful 
within a computational context. Depending on whether the sequence of initial data is specified in terms of the ratio R/M , 



or in terms of the mass M, or of the angular velocity ^l, we can rewrite (|6q) respectively as 



^M— ____ _ M, (70a) 



G X (f7M)2^+2 V R 

= -^Ug) ^(-t) (m) ^''^^ 



, _ 1 [(2^+l)!!]2 f£+iy'+' [GM 



where 

''"^24 (£-1)2^ \J+^) ' ^""'^^ 
and 

X=-j^ pr dr, P=^3- (72) 

In general, the integral in (^2|) needs to be computed numerically, but, in the case of a polytropic equation of state 
P = Kp^ with r 2, it can be computed analytically and, in particular, (Jeffrey 1995) 

''Note that + niQ) < 0. 



12 



Constructing a Mass-Current Radiation-Reaction Force 



1 r.=„. 2.+i._ (2^+1)! 



^2^+3 - 3^2£+i sma;2;-'^dx- ^^^£+1 |E( l^'^' (2£ _ 2A; + 1)! J ' ^^^^ 

Suppose we now impose the condition that the computational timescale , expressed as a multiple N of the stellar 
rotations, be identical to the growth-time t^^^ 

r,^N^^\r,^\ (74) 



Using d?^ ) and (70c), expression ( [7^ ) then becomes a condition on the ratio c?R/ {GM), i.e. 

2/(2£+3) 



c'R 
GM 



0.68 iV^/^ , (75) 



where the numerical coefficient comes from considering £ = 2, F = 2 and 17^ — fi. According to ([Tq), it is always possible 
to rescale the value of the constant c in such a way as to make the growth-time compatible with tEe timescale over which 
the numerical computations can be carried out. Provided we maintain the inequality = \t^^\ 3> ^2"^ or TV ^ 1, 
this rescaling should in no way affect the profiles of the physical parameters during the evolution, but only shorten the 
evolution time over which their growth an d sat uration occur. Alternatively, one can choose M / R to be sufficiently large 
to reduce the growth-time in accord with ( |70bD and then scale the results to stars with more realistic compaction ratios. 
A similar rescaling technique has also been adopted to accelerate the cooling of a hot neutron star and study its collapse 
to a black hole (Baumgarte, Shapiro and Tcukolsky 1996). 

7. CONCLUSIONS 

We have presented a new set of (0 + 3.5) FN hydrodynamical equations in which a 3.5 FN radiation-reaction force due 
to a time-varying mass-current quadrupole moment is considered. Within this system of equations, the hydrodynamics 
is essentially Newtonian except for the inclusion of the relativistic nonconservative effects related to the emission of 
gravitational radiation. 

We have cast this set of equations in a form which is suitable for numerical implementation. In the alternative 3.5 FN 
approach by Blanchet (1993, 1997), the radiation-reaction terms depend on the fifth time derivative of the mass-current 
quadrupole moment Sij . Evaluating such a term accurately within a standard second order numerical scheme could pose 
a problem. Instead, we have chosen a particular gauge in which the radiation-reaction effects depend at most on the 
fourth time derivative of Sij and can therefore be calculated accurately. The additional complication that arises with this 
gauge choice arc three linear, elliptic equations for the components of a vector potential. The solution of such equations is 
no more difficult than the solution of Foisson's equation for the Newtonian gravitational potential and can be performed 
in an identical fashion. 

Simulating the onset and growth of the r-mode instability in rotating neutron stars is highly desirable. A (0 -I- 3.5)FN 
approach may be considerably simpler than a fully general relativistic one and allows one to neglect all conservative 
relativistic effects, which should be perturbative, and focus exclusively on the radiation-reaction effects due to a time- 
varying mass-current quadrupole moment. We have also proposed a suitable rescaling that will make the timescale for the 
onset and saturation of the r-mode instability compatible with any reasonable integration time imposed by computational 
constraints. Work is presently in progress to implement these equations in a numerical code (Ruffert et al. 1999). 

We are grateful to L. Blanchet for his helpful comments and for carefully reading the manuscript. H. Asada would like 
to thank Y. Kojima for useful conversations. This work was supported by NSF Grant AST 96-18524 and NASA Grant 
NAG 5-7152 at Illinois and a JSFS Fellowship to M. Shibata for Research Abroad. M. Shibata also acknowledges the 
kind hospitality of the Department of Physics of the University of Illinois at Urbana-Champaign. 



APPENDIX A 



(iV) 



STRATEGY FOR COMPUTATION OF Sjj 



In this Section we present the analytic integral expressions for the first, the second, and the third time derivative of 
the mass current quadrupole moment Sij. The expressions derived here can then be used, after taking numerical time 

derivatives, to compute S^j^ and (if necessary) S^jK 

The first time derivative of 5"^^ is easily derived from (|l^), after setting Vi = dxi/dt, to yield 



^ki(i J d^x p(^VjfVi + Xj-)^^Xk + 0{e) ; 



S-j^ = <^ki(i I rf'^x p[ Vj)Vi + Xj)'^ ]xk + 0{e) ; (76a) 
= ^ki{i I P {vj)Vi - Xj)di4)) Xk 4- 0(e) . (76b) 



In deriving (76t) we have used the continuity and the Newtonian Euler equation, 



dt '•' dxi 



1 dP d(j) 
p dxi dxi ' 



and exploited the following identity 

efcm(i j rf^X Xj)XkdmP . 

A similar procedure is used for the second time derivative, which can be written as 
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(77) 



(78) 



- / d^x p [x.j)Vkdi(l3 + Xj)Xk [di{dt(f>) + Wm9;,„(/)]] > + 0(e) ; 



(79a) 



efci(i< / d^yiPxk [dj)Vi + diVj)] + / d^x\/ ■ (pv) [xj-^Xkdi 



- / (fyipxk [xj)di{dt(j)) + vidj)(j) + Vj)di<p\ \ + 0(e) 



(79b) 



Note that we have proposed two different expressions for S'^^\ where the second one [i.e. ([79b|)] made use of the following 



identity 



d^x p [v{idj)(l) + VkX(^idj)k(l)\ = j d^x pvkdk [x(^idj)4>] J^d,^^ dk [pvkX(^idj)(t)] - J d^x [x(^,dj)(t>] V • (pv) 

-=/ [x^,d^)4>] V • (pv) , (80) 



in order to eliminate the mixed second partial derivatives of the gravitational potential. 



Expressions (76a)-(79b) apply for a generic equation of state. However, in deriving Slj we will need a time derivative 



of the volume integral of the pressure and a specific equation of state must be specified. In the case of an equation of 
state P = (r — l)pe, we obtain 



c(3) _ , 
- £ki{i 



i^x (r - l)p£|(l - T)xk {dj)Vi + diVj-)) (dmVm) + Vk {dj)Vi + div^) ^ 

+Xk [dj)ai + diaj) - (9j)i'„)(9„u/) - ((9^^;„)((9„UJ))] | 
d^x p {xj)akdi4> + 3vj)Vkdi4> + Xk [aidj)4> + 2aj)di4> + vidj){dt(t>) + viVndj)n(f\ } 
d^x p{[ivj)Xk + 2xj)Vk][di{dtcj)) + fn9/,i0] + Xj)Xk [diduij) + 2vndin{dtcj)) + a„9;„0]} 



d^X Xj)Xkdl{pVmVn)dmn(l> 



0{e) , 



(81a) 



or equally 
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d^xF [d.j)Vi + diVj)] {vk - (r - l)xkdnV„) 

+ f d^xPxk [dj)ai + dlttj) - {dj)Vn){dnVi) - {dlVr,){dnVj))\ 

dl(f) {vj)Xk + Xj)Vk) + Xj)Xk [di{dt(t>) + VndnlC, 



j d^xV • (pv) 



where we have defined 



(fy.Xj)Xkdl(j) [dn {pan) - dn{pVmdmVn)] - / d^XpVk [xj)dl{dt(t)) + {dj)(p) VI + Vj)dl(p] 



(fixpxk 



[2di{dt(j)) + Vndnicj)] + xj^ [di{dtt(p) + v^dniidf 



(81b) 



and used the relation 



d 
dt 



dvj 
'dt 



1 



diP - di(j) , 



J d^xP = J d^x (^^^ = - J d^x P(r - l)djVj 



The partial time derivatives of the gravitational potential (dtcf)) and {dtt't') appearing in (|79| ) and 
following elliptic equations (Nakamura and Oohara 1989) 



(82) 

(83) 
satisfy the 



A{dt^) = -inGdkipvk) 



(84) 



A{dtt^) = ATTG[d.,jipv,v,) + AP + dM(f>)] , 



(85) 



where A denotes the flat spatial Laplacian. While S'^^^' and S^^^ ^ could also be expressed through similar integral ex- 
pressions, this is not useful in general. For a numerical method which is accurate in time and space at the order n, the 
maximum time derivative of Sij which can be calculated accurately is {n + 1). This is because, for a generic S'^"^^'', we 
need n spatial derivatives and n partial time derivatives of (fi. As a result, if one is using a numerical method which is 
second order accurate in space and time, analytic integral expressions are reliable at most up to S^^ . The fourth and 

(3) 

fifth time derivatives need to be obtained by finite differencing of S^^ with increasingly larger truncation errors. This 
consideration is the guideline for our formalism, for which we need only compute S", 



(4) 



APPENDIX B 
CANONICAL FORM OF THE LINEAR METRIC 



Here we obtain the parts of the metric associated with radiation-reaction potential. We follow the notation of Blanchet 
(1993). We first recall that the components of the linearized metric h'^^^ in canonical form and in the harmonic gauge 
condition (Thorne 1980) are given by 



(5 = 0, 

4Gf. (-1)'/ 



-SbL-l (t 

r \ c 



(86) 
(87) 



_8G^ (-1)'/ 



da 



L-2 
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where we have considered only the terms related to mass-current niuhipole moments Sl. Here Ol = di-^di^ ■ ■ • i9i, and 
L — iiZ2 ■ ■ - ii is a compact expression for / indices. Since there is no "tail" term at the linear order, we can derive the 
radiation-reaction metric by taking the half-difference of the retarded and advanced waves (Blanchet 1993) 



1=1 



P - ^2^^ (ITT)! '"^"^"^-^ 



SbL-i{t - rjc) ~ SbL-ijt + r/c) 
2r 



-ab{i 



s'^L-2(t-r/c)-S^^_,(t + r/c) 



2r 



(89) 
(90) 



Expanding the n umer ators of th e right-hand sides with respect to c^^ in order to determine the near zone metric, we 
obtain equations ( |26a| ) and ( p6b| ) as the lowest order of the I ~ 2 mode. The gauge transformation necessary in order to 
set the "new" linear shift = is therefore simply given by 



4G^ (-1)'^ . 

1=2 ^ ' 



2r 



Using (|9l|), the expression of hij in the new gauge is 



(91) 



riij = riij + 2()(i^j) = — 2^ /, , 1 M n; , -, OaL-2 



c- — (? + l)!2/ + l 

1=2 ^ ' 



£ab{i- 



8G 



1=2 



-lyi 



2r 



Sil-\{t~r/c)-Sllt[{t + r/c) 
2r 



(92) 



where Ol is the (symmetric) trace- free part of Ol ■ Equation ( |92D should be compared with the equivalent one obtained in 
Blanchet 's gauge (we recall that we report here only the contributions due to time- varying mass-current multipoles) [cf. 
equation (2.8c) of Blanchet 1997] 



, Blanchct's 
' gauge 



8G^ 

^ 1=2 



{-lyi 21 + 1 



(i + iy. I 



L-1 



Sil]_\{t-r/c)-sill\{t + r/c) 



2r 



(93) 



Expression (p3|), as well as the the second term in (£^), provide no contribution at 3.5 PN order. 

Using (^) and (^2|) at the lowest order of the PN expansion of the I — 2 mode, we obtain equations (E9| 



APPENDIX C 



TIME-SLICE CONDITION AND EQUATION FOR gO 

We here discuss the choice of a time-slice condition and the derivation of the elliptic equation ( ^o|) for ga. We start by 
considering the evolution equation for the 3-metric jij 

^dt%j = -2aK,, + D,l3j + A , (94) 

and the momentum constraint equations 

D,K^^ - DJ< - ^ J, , (95) 

where Di is the covariant derivative with respect to 7^-, Kij is the extrinsic curvature tensor, and K = K^^. The current 
source term on the left hand side of ( p5| ) is defined as = — 7pi/T'^"nc, with — (1/a, —(3^ /a) being the normal to the 
spatial slice. The 3.5 PN expressions of (|9^) and ( ^5| ) are given respectively by 

\dt{jh,^) = -2^K,^ + d^isPj) + ^,{sP^) , (96) 

and 
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where iJi = (Asada, Shibata, and Futamase 1996). 

Taking a further spatial derivative of (Bq) and using the constraints ( 



(97) 



we then obtain 



(98) 



Performing now an infinitesimal gauge transformation yielding [cf. equations (p7|)-(29b)] 



32 G (4) 

7"-y -^^k^klii^jy ■ 

will set to zero the right hand side of equation ( p8| ) [we recall that di{7hij) 
a result, we choose as slice condition 



(99) 
(100) 

0] and thus require to be a constant. As 



K = 



(101) 



at all times. Condition ( |101[ ) is known as the maximal time-slice condition (Arnowitt, Dcscr and Misner 1962; Smarr and 
York 1978; Schafer 1983; Blanchet, Damour and Schafer 1990). As a consequence of (101), the evolution equation of K 
is given by 



1 



dtK 



where = T^.^nf^n" and S = T^.^-f'"' . 



The metric coefficient ga is then obtained as the solution of the 3.5 FN expression of the elliptic equation (102) 

AttG 



-a{p^+S) + aK.jK'^ 



(102) 



(103) 



After discarding all the 3.5 P N ter ms except those arising from the mass-current quadrupole, the 3.5 FN expression of 
the left-hand side of equation (103) is written as [cf. equation (g)] 



while the (full) right-hand side of equation ( |103| ) is rewritten as 



47rG , 
— ^a(p^ 



S) 



AirGa 



1 1 I P 

1 + - £ + - 
c \ P 



1 + —f^W.Wj + —P 



(104) 



(105) 



It is easy to estimate that, at the 3.5 FN level, the contributions to expression ( |1G5| ) coming from the mass current 
quadrupole moment are at most 0{c~^^) and that the slicing condition for ga is therefore given by {M. This follows 
from the fact that the contribution in p, e, and P is 0(1), and is 0(c~^) for •y'^^WiWj . As a consequenceTthe contribution 
of the mass current quadrupole moment from the terms in the cur ly b rackets of ( |1G5| ) is at most of 0{c~^). Similar 
considerations apply also for the last term in the right-hand side of (|105| ) where Kij = 0{c^^), while the contribution of 
the mass-current quadrupole moment in Kij is 0(c~^) [cf. equation (M)]. As a result, the mass-current contribution in 
aKijK^i is at most 0(c-"). 
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